This study uses Nihommatsu city and Mimanisoma City as case studies. The two cities are 70km and 30km from Fukushima Daichi Plant respectively.
library("leaflet")
library(readr)
library(dplyr)
library(RColorBrewer)
library(Hmisc)
NB:Important Notes from the Nuclear Regulation Authority on JAEA website
Purposes of this vehicle survey were: 1. To ascertain the tendency and cause of time change of air dose rates by comparing past vehicle survey data and survey meter data at the height of 1 m above ground as well as “walk survey” data, and 2. To contribute to the establishment of radioactive substances distribution prediction model. MEXT evaluated the decrease in the air dose rates caused by the decay of cesium during the survey period and it was around 1%, which was smaller than the errors of measuring instruments source
niho <- read_csv("jun_2011_fukushima.csv")
niho2013 <- read_csv("nov_2013_fukushima.csv")
# View(niho2013)
# names(niho)
dim(niho)
## [1] 45273 18
# head(niho)
class(niho)
## [1] "tbl_df" "tbl" "data.frame"
names(niho) <- c("gridcode","pref","city","gridCenterNorthlat","gridCenterEastlng","gridCenterNorthlatDec",
"gridCenterEastlngDec","daichi_distance","no_samples","AvgAirDoseRate",
"NE_nLat","NE_eLong","NW_nLat","NW_eLong",
"SW_nLat","SW_eLong","SE_nLat","SE_eLong")
names(niho2013) <- c("gridcode","pref","city","gridCenterNorthlat","gridCenterEastlng","gridCenterNorthlatDec",
"gridCenterEastlngDec","daichi_distance","no_samples","AvgAirDoseRate",
"NE_nLat","NE_eLong","NW_nLat","NW_eLong",
"SW_nLat","SW_eLong","SE_nLat","SE_eLong")
#Strip Nihommatsu city,
niho$city[niho$city == "Nihommatsu city"] <- "nihommatsu"
niho2013$city[niho2013$city == "Nihommatsu city"] <- "nihommatsu"
nihom <- subset(niho, city == "nihommatsu")
nihom2013 <- subset(niho2013, city == "nihommatsu")
dim(nihom)
## [1] 2944 18
dim(nihom2013)
## [1] 11758 18
# View(nihom2013)
class(nihom)
## [1] "tbl_df" "data.frame"
niho_q <- nihom %>%
mutate(dose_quants = cut2(nihom$AvgAirDoseRate,cuts=c(0.1,0.5,1.0,1.5,2.0,2.5,3.0),levels.mean=TRUE))
# View(niho_q)
niho_q <- na.omit(niho_q)
# write_csv(niho_q, path = "niho_q.csv")
nihom2013_q <- nihom2013 %>%
mutate(dose_quants = cut2(nihom2013$AvgAirDoseRate,cuts=seq(0.06,1.6,0.25),levels.mean=TRUE))
summary(nihom2013$AvgAirDoseRate)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0630 0.2800 0.3700 0.4009 0.4900 1.4000
summary(niho_q$AvgAirDoseRate)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.190 0.710 1.000 1.029 1.300 2.500
iro <- colorFactor(
palette = "YlOrRd",
domain = niho_q$dose_quants
)
iro2013 <- colorFactor(
palette = "YlOrRd",
domain = nihom2013_q$dose_quants
)
# Link of Daichi
fukulink <- paste(sep = "<br/>",
"<br><a href='http://www.tepco.co.jp/en/decommision/index-e.html'>Fukushima Daichi</a></b>",
"Source of radiations"
)
niho_plot <- leaflet() %>%
addTiles()%>%
addRectangles(data = niho_q,lng1 = ~SW_eLong, lat1 = ~SW_nLat,lng2 = ~NE_eLong, lat2 = ~NE_nLat,
color = ~iro(niho_q$dose_quants))%>%
addLegend("bottomright", pal = iro, values = niho_q$dose_quants,
title = "AvgAirDoseRates",
labFormat = labelFormat(prefix = "µSv/h "),
opacity = 1)%>%
addPopups(lat = 37.4211, lng = 141.0328,popup = fukulink,
options = popupOptions(closeButton = TRUE))
niho_plot
niho2013_plot <- leaflet() %>%
addTiles()%>%
addRectangles(data = nihom2013_q,lng1 = ~SW_eLong, lat1 = ~SW_nLat,lng2 = ~NE_eLong, lat2 = ~NE_nLat,
color = ~iro2013(nihom2013_q$dose_quants))%>%
addLegend("bottomright", pal = iro2013, values = nihom2013_q$dose_quants,
title = "AvgAirDoseRates",
labFormat = labelFormat(prefix = "µSv/h "),
opacity = 1)%>%
addPopups(lat = 37.4211, lng = 141.0328,popup = fukulink,
options = popupOptions(closeButton = TRUE))
niho2013_plot
ggplot(niho_q, aes(daichi_distance,AvgAirDoseRate)) +
geom_point() +
geom_smooth(se = FALSE)+
ggtitle("AvgAirDose against Distance to Daichi Plant")
ggplot(data = niho_q) +
geom_bar(mapping = aes(x = daichi_distance, fill = dose_quants), width = 1)+
ggtitle("AvgAirDose Measured Counts against Daichi Distance")
mina <- read_csv("niho.csv")
dim(mina)
## [1] 45273 18
# View(mina)
#change Minamisoma city to minamisoma
mina$city[mina$city == "Minamisoma city"] <- "minamisoma"
#filter nihomatsu observations only
mina <- subset(mina, city == "minamisoma")
dim(mina)
## [1] 1985 18
summary(mina$AvgAirDoseRate)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.210 0.520 0.780 1.771 1.800 18.000
#plot(mina$AvgAirDoseRate)
mina_q <- mina %>%
mutate(mina_quants = cut2(mina$AvgAirDoseRate,cuts=seq(0.2,20,2.2),levels.mean=TRUE))
## View(mina_q)
mina_q <- na.omit(mina_q)
write_csv(mina_q, path = "mina_q.csv")
mina_iro <- colorFactor(
palette = "Set1",
domain = mina_q$mina_quants
)
mina_plot <- leaflet() %>%
addTiles()%>%
addRectangles(data = mina_q,lng1 = ~SW_eLong, lat1 = ~SW_nLat,lng2 = ~NE_eLong, lat2 = ~NE_nLat,
color = ~mina_iro(mina_q$mina_quants))%>%
addLegend("bottomright", pal = mina_iro, values = mina_q$mina_quants,
title = "AvgAirDoseRates",
labFormat = labelFormat(prefix = "µSv/h "),
opacity = 1)%>%
addPopups(lat = 37.4211, lng = 141.0328,popup = fukulink,
options = popupOptions(closeButton = TRUE))
mina_plot
ggplot(mina_q, aes(daichi_distance,AvgAirDoseRate)) +
geom_point() +
geom_smooth(se = FALSE)+
ggtitle("AvgAirDose against Distance to Daichi Plant")
ggplot(data = mina_q) +
geom_bar(mapping = aes(x = daichi_distance, fill = mina_quants), width = 1)+
ggtitle("Counts of Measuring Locations per Air Dose Rate")
library(ggplot2)
niho_q <- read_csv("niho_q.csv")
niho2013_q <- read_csv("niho2013.csv")
names(niho_q) <- c("gridcode","pref","city","gridCenterNorthlat","gridCenterEastlng","gridCenterNorthlatDec",
"gridCenterEastlngDec","daichi_distance","no_samples","AvgAirDoseRate2011",
"NE_nLat","NE_eLong","NW_nLat","NW_eLong",
"SW_nLat","SW_eLong","SE_nLat","SE_eLong")
names(niho2013_q) <- c("gridcode","pref","city","gridCenterNorthlat","gridCenterEastlng","gridCenterNorthlatDec",
"gridCenterEastlngDec","daichi_distance","no_samples","AvgAirDoseRate2013",
"NE_nLat","NE_eLong","NW_nLat","NW_eLong",
"SW_nLat","SW_eLong","SE_nLat","SE_eLong")
niho11_13 <- merge(niho_q, niho2013_q, by.x = "gridcode", by.y = "gridcode", all = TRUE)
# View(niho11_13)
# Check if the merged columns are real identical
niho11_13 <- na.omit(niho11_13)
identical(niho11_13$gridCenterNorthlat.x,niho11_13$gridCenterNorthlat.y)
## [1] TRUE
identical(niho11_13$daichi_distance.x,niho11_13$daichi_distance.y)
## [1] TRUE
identical(niho11_13$no_samples.x,niho11_13$no_samples.y) #2011 and 2013 samples differ,but lets ignore that.
## [1] FALSE
niho11_13 <- select(niho11_13,gridcode,pref.x,city.x,gridCenterNorthlat.x,gridCenterEastlng.x,gridCenterNorthlatDec.x,
gridCenterEastlngDec.x,daichi_distance.x,no_samples.x,AvgAirDoseRate2011,
NE_nLat.x,NE_eLong.x,NW_nLat.x,NW_eLong.x,no_samples.y,AvgAirDoseRate2013)
# create new data set (niho11_13)
write_csv(niho11_13,path = "niho11_13.csv")
niho11_13 <- read_csv("niho11_13.csv")
#View(niho11_13)
#View the change
plot(x=niho11_13$AvgAirDoseRate2011,type="l",col="red")
lines(niho11_13$AvgAirDoseRate2013,col="green")
#Frequency
hist(niho11_13$AvgAirDoseRate2011,col="green")
hist(niho11_13$AvgAirDoseRate2013,col="red",add=TRUE)
continuous interactions (two continuous variables) with model \[Y{i} = \beta_{0} + \beta_{1} X_{1i} + \beta_{2} X_{2i} + \beta_{3} X_{1i} X_{2i} + \epsilon_{i}\]
AvgAirDoseRate2013 = {0} + {1} AvgAirDoseRate2011 + {2} daichi_distance.x + {3}(AvgAirDoseRate2011)(daichi_distance) + _{i}$$
fit1 <- lm(AvgAirDoseRate2013~daichi_distance.x, data = niho11_13)
summary(fit1)
##
## Call:
## lm(formula = AvgAirDoseRate2013 ~ daichi_distance.x, data = niho11_13)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.27317 -0.10516 -0.02457 0.07873 0.89017
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5844396 0.0200218 29.19 <2e-16 ***
## daichi_distance.x -0.0049126 0.0003721 -13.20 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1373 on 2805 degrees of freedom
## Multiple R-squared: 0.0585, Adjusted R-squared: 0.05817
## F-statistic: 174.3 on 1 and 2805 DF, p-value: < 2.2e-16
confint(fit1)
## 2.5 % 97.5 %
## (Intercept) 0.545180562 0.623698590
## daichi_distance.x -0.005642192 -0.004182925
ggplot(niho11_13, aes(x = AvgAirDoseRate2011, y = AvgAirDoseRate2013)) +
geom_point() +
stat_smooth(method = "lm", col = "red")+
ggtitle("Correlationg between Air Dose Rate of 2011 vs of 2013")
End